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ABSTRACT 


Long term numerical integrations were performed on a 
two-dimensional model utilizing the barotropic vorticity 
equation including a space-and time-dependent forcing func- 
tion and a laminar-type viscosity term. Three different 
forms of the laminar type viscosity or dissipation term 
were used. A comparison was made of the effects these 
different dissipation terms had on the time-dependence of 
total kinetic energy, total enstrophy, and total gradient 
of vorticity squared, the spectral distribution of kinetic 
energy and enstrophy, and the stream function field. - A 
difference of behavior of all the above variables was found 


for each form of dissipation. 
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I. a INERODUCTITORK 


With modern computers it is possible to simulate turbu- 
lent motion in three dimensions., As a consequence of theo- 
retical studies as well as experimental results, the 
turbulent processes at work in the three-dimensional case 
are reasonably well understood. In recent years, however, 
Studies have also been conducted on two-dimensional turbu- 
lence as an analogue to the large scale motions of the atmos- 
phere. The behavior of the two-dimensional turbulence 
parameters is just beginning to be discovered. 

The principal distinguishing feature of two-dimensional 
turbulence is the tendency of its energy-containing eddies 
to grow indefinitely in scale. In any finite difference 
model of two-dimensional turbulent flow, therefore, the 
energy transfer between large scales will be computed more 
or less correctly. At the scale of the grid interval, how- 
ever, aliasing errors can transfer energy erroneously to all 
of the larger scales destroying the sense of the calculation. 
This problem can be alleviated by the removing of energy 
from scales near the grid scale. The energy thus "dissipated" 
is that which would cascade to smaller scales in real flow. 

A common procedure in numerical modeling is to para- 
meterize the sub-grid scale dissipation in the form of a 


diffusion process with an assigned coefficient of "eddy 





diffusion", A, as in V °> (AVS), where cea open 
being dissipated. 

The objective of the present work is to compare the 
behavior of certain properties generated by different forms 


of the eddy diffusion term in two-dimensional turbulent flow. 
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Il. THE, MATHEMATICAL MODHE 


The principal equation of the model is the barotropic 
vorticity equðťion including a space- and time-dependent Miorre 


ing function and a laminar-type viscosity, written as 
3L iu 
XEM ME M we) @ 


The wind velocity is associated with the stream function 

by the following equation: 
we gr + SÈ (2) 
The vorticity is the lapłacian of the stream function, 
S ES (3) 

In all of the numerical experiments undertaken in this 
study, the regime of integration is taken to be a rectangle 
kith cyclic boundary condition in the x direction and free 
slip boundary conditions on the northern and southern bound- 
aries. (See Appendix A). 

The numerical integration of Equation (1) and (3) was 
performed by first replacing them with a pair of finite dif- 
ference equations to be solved for the variables prescribed 
On a cartesian grid. (See Appendix B). 

The numerical space-differencing method used on the non- 
linear terms of Equation (1) is the second-order scheme 
developed by Arakawa (1966) which conserves kinetic energy 
and enstrophy within the nonlinear terms. The laplacians 


.in Equations (1) and (3) were approximated by the methods 


1l 





discussed in Appendix B. Time integration was performed by 
the explicit centered leapfrog scheme, Bene for the first 
step and restarts every 50 time steps by the iterative Euler- 
backward scheme. (See Appendix D). 

The forcing function S used in these experiments was a 
crude approximation for the baroclinic driving forces in the 
atmosphere. S was negatively related to the stream function 
and was scaled to the magnitude of the rate of change of 
vorticity. (See Appendix A). The source term was set to 
zero at mid-time in each experiment conducted to study the 
decay properties of each form of dissipation. 

The initial stream field was defined to have components 
in wave numbers 2, 4, and 6 and of such magnitude as to 
simulate the 500-mb patterns. (See Appendix A). 

The different forms of the dissipation are discussed in 


the following chapter. 
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Lil. FORMSTFOF DISSIFATLON 


In Equation (1), the last term on the right hand side 
represents a laminar-type viscosity or dissipation and is 
written as V * (AVz). The experiments assumed three differ- 
ent forms of A which represents the coefficient of "eddy 
diffusion." 

In the first form, A was taken to be the constant .5 x 
de e Se This value was used to produce an equilibrium 
solution of a one-wave stream function field within the time 
frame of these experiments. 

The second form of A was taken to be a function of the 
wind field as formulated by Smagorinsky (1963) for three- 
dimensional turbulence: 


Z 
A. = Y, |D,| (Ax) (4) 


muere D. is the deformation on the grid Ax and Yi is a con- 
stant. 


The deformation is defined as: 


[D cay ae De 
where 
- ou _ Iv 
PES X dY 
and 
DINO NI AO 
Ds ðX i oy 


The third and final form of A was that proposed by Leith 


(1969) for two-dimensional turbulence, namely, 


I3 





Ay = Yo M4 (Ax)? (5) 


Here Y» is a constant and | Vz | is the scalar magnitude 
of the cradient of vorticity: 
The finite difference forms of the above equations may 


be found in Appendix B. 
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IV. DISCUSSION OF RESUCES 


Seven experiments were conducted for comparison and 
evaluation of the different forms of dissipation as follows: 
(1) . Equation (1) with the source and dissipation terms 
removed. 
(2) A constant coefficient for A in Equation (1). 
(3) Smagorinsky's form of A (Equation 4) with Y ks 


(4) Smagorinsky's form of A with qe SE 


(5) Smagorinsky's form of A with Yi 13% 

(6) Leith's form of A (Equation 5) with one $94 

(7) Leith's form of A with E 22" 

In each experiment, the time-dependence of total kinetic 
energy, total enstrophy, and the total magnitude of the 
gradient of vorticity squared were calculated, along with 
the spectral distribution of energy and enstrophy. The 
effects of each ioe ees. on these measurements and 
on typical wave disturbances were of prime interest. 

Comparing first the time-dependence of kinetic energy 
(Figures 2, 3, 4), it was noted that except for the zero dis- 
sipation, zero source case there was a tremendous decrease of 
kinetic energy initially. The order in which this dissipa- 
tion changed into regeneration was, by experiment numbers, 
7,2,3,5,4,6. In general the Leith form Nalted= eters inter 
rapid dissipation of kinetic energy earlier. 


The energy in the model at mid-time was, by decreasing 


order of magnitude, 1,5,7,2,4,3,6.  —<These valtes weee 
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obviously dependent upon tlie size ofthe dissipation tern: 
However, the energy present at end-time was not in the same 
order as above, 1,5,7,4,3,2,6. This indicated the A = con- 
stant form had a greater rate of dissipation than the other 
forms after the source term was taken out. An interesting 


point noted was that for both A, and A 


L e? the smaller Y terms 


produced the greater rate of dissipation. This was due to 


the dependence of both A, and As on the strength of the flow 


L 
and with the smaller y the flow was stronger as will be 
discussed later. 

Looking next at the time-dependence of enstrophy (Figures 
5,6,7), in each of the experiments it was easily seen that 
enstrophy was far more conservative than energy. However, 
the enstrophy dissipated rapidly initially as did the kinetic 
energy. In the Smagorinsky and Leith forms of A there was a 
tendency for the enstrophy to stabilize faster for larger 
values of Y. Leith's form caused a slight growth in enstrophy 
after the initial drop but once the source was cut off at 
mid-time enstrophy dissipated slightly. In Smagorinsky's 
form, enstrophy was more nearly conserved while the source 
was in Equation (1), but enstrophy had a tendency to be dissi- 
pated when the source was turned off. A constant A did not 
conserve enstrophy. 

The final time-dependent term to be considered was the 
total gradient of vorticity squared (Figures 8,9,10). This 
term was even more conservative than the two terms considered 


before. There was a tendency for equilibrium with respect 
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to |vz|?to be reached faster the greater the Y terns. 

The fastest equilibrium solution was from experiment 6 with 
experiment 3 following. The experiments in which the total 
gradient of vorticity squared was most conserved were both 
Leith forms and the Smagorinsky form with Y=.5, while there 
were a number of small fluctuations in 2. The deviations 
noted in the gradient of vorticity squared terms seemed to be 
more like step functions than smooth curves after the initial 
level-off of values in all experiments other than 1. 

The next measurements to be considered were the spectral 
distributions' of kinetic energy and enstrophy (Figures 11-18). 
The scalar energy spectrum and the scalar enstrophy spectrum 
were plotted on logarithmic scales at the initial time, 1512 
time steps, and 3024 time steps. At the 1512th time step the 
Source term was taken out of Equation (1). The 3024th time 
step was the final time step. The values of kinetic energy 
and enstrophy used in the spectral analysis were taken off 
the middle grid row, K=8. Since there were 25 grid points 
east-west, it was possible to obtain 12 wavenumbers. Con- 
ponents with wavenumbers greater than 6 are subject to 
possible distortion from aliasing errors. 

In the initial spectra of kinetic energy and enstrophy, 
the energy and enstrophy were mainly contained in the even 
Wavenumbers, which was due to the form of the initial stream 
field, and evenly distributed in these wavenumbers. 

The zero dissipation, zero source experiment spectra had 


no discernible pattern other than a tendency for both energy 
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and enstrophy to be spread fairly evenly in all wavenumbers. 
Experiment 2, where A was a constant, concentrated energy 
and enstrophy in the lower even wavenumbers without any dis- 
tinctive pattern., 

In experiment numbers 3, 4, and 5, there was a tendency 
for a -3 power slope to become apparent in the plots of the 
spectral variation., It was noted in the Ac form, as Y became 
smaller the spectral distribution tended to "bend" away from 
the -3 power slope. The dropoff at the large wavenumbers end 
of the spectrum also became more pronounced as y became smal- 
ler. Obviously, the amount of energy and enstrophy present 
in corresponding wavenumbers was increased with decreasing Y. 

Experiment numbers 6 and 7 had a sawtooth appearance with 
the even numbered waves generally containing more energy than 
the odd numbered waves. The envelope of the spectral distri- 
bution did have a tendency toward the -3 power slope, with the 
strongest tendency one na n mid-time.  Enstrophy tended 
to approach the -3 power slope more nearly than kinetic energy. 

For a comparison of the effects of the different dissipa- 
tion forms on the wave disturbances, the stream function field 
was analyzed at the initial time and after time steps 756, 
1512, 2268, 3024 or approximately every 5 days (Figures 19-47). 
All experiments were started with the same initial y field. 
The zero dissipation, zero source experiment may be referenced 
in the figures., 

After time step 756, the stream function field was begin- 


ning to become quite zonal in character. It appeared to have 
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primarily long wave components that were dominated by wave- 

number 2. The strongest flow pattern was in the experiments 
with the smallest y and in the A = constant case. Low areas 
were situated in the zonal band around the middle grid row, 

K=8, High areas were along the north-south boundaries. 

At time step 1512, mid-time for all experiments, the flow 
pattern was very zonal. The low areas again were along the 
middle grid rows thereby causing easterly flow in the north 
and westerly flow in the south. The experiments with the 
meet zonal flow were 7, 5 and 2. Leith's form of dissipation 
with y=.5 produced the least zonal flow with a wave form more 
obvious than the other forms of dissipation, 

After time step 2268, the flow is still predominately 
zonal; however, it had become weaker and was starting to 
deform in all experiments, the most notable cellular circula- 


pron being in A, with y=.5. 


L 
At the end-time of the experiments, the flow had con- 
siderably weakened.  Leith's form of dissipation, with y=.2 
and the A=constant form stream function patterns were still 
very zonal in nature. Smagorinsky's forms deviated from the 
zonal flow more as the Y term decreased. The Smagorinsky 
form with Y17:13 had cellular circulation corresponding to 


about wavenumber 2. Leith's form with y=.5 had the weakest 


flow pattern and the least zonal. 
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V. CONCLUSIONS 


Leith"s form of “dissipation was the fastest (ecm NC 
the dissipation forms. All of the time-dependent curves 
reached their minimum value faster with Ar After the dis- 
sipative effect had leveled m the source term had become 
dominant and the Ar forms had more of an opportunity to 


grow, with the smallest y coefficients growing most rapidly. 


The spectral curves showed that A, had a tendency to form an 


L 
envelope which approximately fit the -3 power slope. The 
variation in the constant y in Leith's form of dissipation 
had the greatest influence on the stream field. 
Smagorinsky's form for dissipation was not as fast act- 
ing as Leith's form. The Smagorinsky's form was also not 


as sensitive to the constant coefficient as was A lvo 


L? 
all measurements made with different y's indicated well be- 
haved curves that related almost one to one with each other. 
For the Ag form, the spectral distributions between Wave- 
numbers 1 to 6 fit the -3 power slope fairly well, with 
TS, 2 being the best fit. The stream function field for Ag 
was generally the least zonal of the forms of dissipation 
considered, with the strength of the flow inversely dependent 
upon the constant coefficients. The more zonal the flow, the 
Stronger the gradients of the stream function appeared. 

The constant form of A had the slowest reaction time in 


the time-dependence measurements as long as the source term 


was still in Equation (1). However, Ac produced the greatest 


20 





variations when the source term was removed in the time- 
dependence measurements. In the spectral measurements, the 
Ao form tended to have the -3 power slope in the even 
numbered waves only. The stream function field produced by 
Ab moved toward zonal flow at about the same rates as the 
other forms; however, once the zonal flow was established, 
A form kept the zonal flow more than the other forms. 

Due to the simplicity of the model used in this study 
and the use of synthetic initial data, it is impossible to 
conclude which form of dissipation is the best in meteoro- 
logical terms. However, this author believes that the 
differences in the experiments performed are great enough 
to justify further studies of this nature. 

Additional experiments that appear desirable are the 
following: 

First the expanding of the grid to investigate more 
wavenumbers. 

Secondly, the addition of a linear drag term on the 
right hand side, namely, EC, to simulate chef surfacen uu 
tion forces at work in the atmosphere. 

Thirdly, and probably the most important, changing the 
source term to better simulate the baroclinic processes of 
the real atmosphere, with the source term having an effective 
wavenumber between 7 and 8. 

It is expected that these variations to this model would 


produce a more realistic spectral distribution, and Ines 


ance with two-dimensional turbulence theory. Hopefully, 
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the new formulation would produce a stream function field 


that was meteorologically meaningful. 
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APPENDIX A 


COMPLETE DESCRIPTION OF THE MODEL 


15 + + + + + + + + + + + ++ + ++ ++ + +++ + G Roc 


13 + + + + + + + + + + + + + + ++ + ++ +++ ++ + 
12 + + + + + + + + + + + + + + + + + + + + + + + + + 
ll + + +. ++ + + +++ ++ +++ +++ + +++ +++ 
10 + + + + + + + + + + + + + + + + + + + + + + + + + 
9 +++ +++ ++ ++ +++ + + + + +++ + + + + + 
8 + + + + + + + + + ++ + +++ +++ +++ +++ 4 
7+ tet ttt ttt ttt tte tte L K L ttt + + 
6 + + + + + + +++ ++ ++ ++ + +++ ++ ++ + + 
5 E +t+t+tatteitteit tetti ttt ++ +4 + 
4 + + +++ +++ +++ ++ +++ ++ ++ +++ + + 
3 +++t+ttettittt+t + ttt + +++ +++ + + + 
do PT 0 mp 
1+ +++ ++ ++ ++ +++ ++ ++ +++ ++ + + + 

pesi 3 4 5 6.7.8 91011 12 13 14 15 16.17 18.19 20 21) 22523 24525 


Figure 1. Grid for the model. 


North-south direction has 15 grid points indexed with K. East-west 
direction has 25 grid points indexed with J. Distance between grid 

: O : 
points was 2.5 latitude or 150 nautical miles (Ax=Ay=2.5°). The grid 


is imposed on a constant P plane located at 45°N (K-8). 


BOUNDARY CONDITIONS 
The southern boundary was placed at K=2 and the northern at K=14. 
The stream function was set equal to zero along these boundaries, i.e. 
V(3,2) 2» ¥(J,14) = 0 
To insure free slip on the northern and southern boundarics 
¥(J,1) = = $(J,3) 
and ¥(J,15)= - (5,13) 


was imposed on the model. 
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Periodicity was assumed in the east-west direction and achieved 
by: 
V(1,K) = Y(26,K) 


and ¥(25,K)= Y(0,K) 


INITIAL CONDITIONS 


The initial stream function was defined as: 


V(J,K) 7 - C sin CB sin Ee - C sin EE) sin = 





Ba Sin E 


L 
where 
x = (J-l) Ax Dea de ad 
Y = (K-3) Ay K = 3,4,...,16 
W = 12 Ay (width between north-south boundaries) 
L = 2W dx (length of grid periodicity) 
C = 1.0 x 10° 


The source term was related to the stream function by 
— | 
S So y /C 


where So = 1.0 x 10719. 
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APPENDIX B 


FINITE DIFFERENCE EQUATIONS 


Equations (1), (2), and (3) are solved numerically by introducing 
finite-differences in x, y, and t. Sequential over-relaxation was 
used to solve for the Y field and is discussed in Appendix C. Time 
differencing is discussed in Appendix D. 

The vorticity, or the Laplacian of Y, was approximated by the 
second-order expression using the center and four immediately conti- 
guous points. 

c <V 2 [E OFLEH (O-L,- V(O,EH)* ((,K-2) -46 (,8)] / Us) ^9 (so 
The numerical space-differencing method used on the non-linear 
term on the left hand side of Equation (1) is the second-order scheme 

developed by Arakawa (1966), which conserves kinetic energy and 


enstrophy within the nonlinear term. 


J(Y,6) ved (Co 4) 
J(C,¥) E (Licor a = CO) ene E V (J,K-1)] 


“IC KH) - CO,K-]. EEG). - ALO] /40* 

* [cone [Y (J,K+1) -9 (3+1, K) ]- IA A T 
A Rene S POL) 
JA (Ic) * 
A C(J-1,K) [Y(J-1,K#1)-Y(J-1,K-1)] 
- C(J,K+1) [Y(J+1,K+1)-Y(J-1,K+1)]+ CQ, K-2 I Q8, k- 1-6 -1, 01) 
/4 (dx)? 

The linear 8 term in Equation (1) was approximated by the first order 


difference equation 


B a = P (Iron E !O-1,91/25s] 
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The above procedure was also used to find the finite difference form 
of the wind velocity components. 

For the final term in Equation (1), an imaginary staggered grid 
had to be imposed. This half grid was needed because the finite dif- 
ferencing procedure limited the values obtainable. The form of the 
viscosity term was 
V* (AVC) = Irun: RER) | [ACI + ACK) 1/24 

- [C(J,K)- C(J-1,K)]/4x | tao, + A(J+1,K)]/2 |) 
+ ME CODIA] | TAGE) AI] 
- | [C(J,K)- C(J,K-1)]/Ay | aus dr vis) 
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APPENDIX C 


RELAXATION TECHNIQUE 


The equation to be solved for the OY/0t field from Equations (1) 


and (3) is 
CR oe Gy) (C-1) 
where 
Py) = IGM) - P Ù+ s +y (QVO, (6-2) 


subject to the boundary conditions, dY/dt = 0 on the north-south 
boundaries and periodic in the east-west direction. 

In the numerical experiments, F(J,K) was computable on all interior 
grid points and the boundary conditions were known. Therefore, values 


of OJ/Ot(J,K) were determinable from the Poisson equation 


vot (J,K) = F,(J,K) (G=3) 
where 
F, (JK) = finite difference form of F(J,K) (Ax) * (C-4) 


The first step in the iterative technique of sequential over-relaxa- 
tion for the solution of Equation (C-3) was to make a first guess of 
values of dV/dt for every grid point to satisfy Equation (C-3). In 
these experiments this first guess after one step was zero everywhere; 
in subsequent steps the value from the previous time step was used for 
the guess field. 

For every estimate of the solution, the residual is defined as 
follows: 


(N) 
ig a Qu (J,K) = F, (3,0 


where the superscript denotes the iteration. 
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A new estimate was made according to the formula: 
N+1) (N) 
at) QU 1+a) _N 2 
t (J,K) = t (J,K) + È R (âx) 


where Q is the overrelaxation coefficient and for these experiments 


equaled .65. 
The iterations were continued until ip eK es © fOr val ede 


S -1 
m 


where € was taken to be 1.0 x 10 sec at which time the solution 


was considered to be sufficiently accurate. 
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APPENDIX D 


TIME STEP SCHEMES 


In general, centered time differences (leapfrog) were used for all 
quantities; however, the first step was a Euler-backward time step. 
The latter was also introduced every 50 time steps to selectively damp 
high frequency waves and to prevent the separation of solutions at even 
and odd time steps caused by the leapfrog time scheme. 


LEAPFROG SCHEME 


E (t-1) 


(t) 
(J,K) + 20t n (J,K) 


(J,K) = Y 
where superscripts denote the time level. 
EULER-BACKWARD 
The Euler-backward is composed of two steps: 
1. (forward step) EE = y C ouo E O eaw 


ay CETL) 


2. (simulated backward) Ta = y uos Ate (J,K) 


Time steps were lO minutes and a total of 3024 steps (about 3 weeks) 


were taken in each experiment. 
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